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Chapter  1.  Int  roduct ion 


Recently,  a  series  of  papers  and  reports  (Yue  and  Met,  1980;  Kirby,  1983, 
1984;  Kirby  and  Dalrymple,  1983a, b,  1984;  Liu  and  Tsay,  1984)  have  detailed 
the  derivation  of  a  parabolic  model  for  the  forward  scattered  component  of  a 
train  of  time-periodic,  weakly  nonlinear  Stokes  waves.  In  its  most  general 
form,  the  model  is  applicable  to  the  study  of  wave  propagation  through  regions 
with  both  slowly  varying  depth  and  ambient  currents.  For  the  case  without 
currents,  good  agreement  between  experiment  and  theory  has  been  demonstrated 
(Kirby  and  Dalrymple,  1984  and  Liu  and  Tsay,  1984)  for  cases  where  Stokes 
theory  is  strictly  valid. 

The  development  of  a  wave  model  designed  for  practical  application  to 
realistic  modelling  problems  poses  a  number  of  physical  problems  which  do  not 
fit  wf  hin  the  theoretical  context  of  the  Stokes  wave  model.  Several  of  these 
problems  are  discussed  in  the  present  report,  and  methods  for  including  their 
effects  in  the  combined  refraction-diffraction  model  are  proposed. 

Tn  Chapter  2,  the  problem  of  matching  dispersion  relations  between  Stokes 
and  shallow  water  waves  is  discussed,  and  a  formulation  is  proposed  which 
smoothly  matches  the  dispersion  relation  for  Stokes  waves  to  an  approximate 
dispersion  relation  due  to  Hedges  (1976)  for  shallow  water.  The  proposed 
matching  alleviates  the  need  for  extending  approximate  calculations  into 
deeper  water  when  their  use  in  shallower  water  is  dictated.  Further,  the 
nrovision  of  a  continuously-varying  formula  for  all  water  depths  insures  a 
smooth  variation  in  both  phase  and  group  velocities  from  deep  to  shallow 
water,  eliminating  the  regions  of  discontinuity  between  Stokes  theory  and  the 
full  cnoidal  wave  theory. 


« 
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In  Chapter  3,  a  formulation  for  waves  propagating  at  not-so-small  angles 
to  the  preferred  direction  and  in  the  presence  of  an  ambient  current  is  given 
in  full.  The  present  derivation  follows  the  scheme  indicated  by  Kirby  (1983), 
but  the  material  presented  here  was  not  pursued  to  the  point  of  completeness 
in  the  former  report. 

Final lv.  Chapter  4  presents  an  application  of  the  breaking  wave  model 
described  bv  Kirbv  (1983)  and  Dally,  Dean  and  Dalrymple  (1984)  to  the  case  of 
a  surf  /.one  aroun  !  an  island.  Since  this  chapter  is  in  the  form  of  a  full 
manuscript,  discussion  of  the  physical  problem  is  left  to  the  chapter 
i ntroduct ion. 


Chapter  2.  A  Modification  of  the  Nonlinear  Dispersion 
Relation  for  Shallow  Water  Waves 

A  central  problem  in  present  efforts  at  modelling  waves  in  coastal 
regions  consists  of  the  mismatch  in  wave  properties  as  waves  propagate  from  an 
intermediate  depth  region,  in  which  Stokes  theory  is  valid,  into  the  nearshore 
zone,  where  waves  ire  better  described  by  the  Boussinesq  equations.  Several 
theories  tor  large  amplitude  waves  in  arbitrary  depth  (Dean,  1965;  Rienecker 
and  Fenton,  1981,  for  example)  provide  a  bridge  across  the  division  between 
Stokes  theory  and  cnoidal  theory;  these  theories  are  computationally  intensive 
and  are  typically  used  to  accurately  describe  the  properties  of  a  single 
wave.  For  the  computation  of  waves  propagating  through  large  areas  of 
variable  depth,  the  use  of  the  theories  for  waves  of  small  amplitude  Is 
computationally  more  feasible  at  present;  formulations  of  this  type  are  the 
subject  of  the  present  report. 

2.1.  Review  of  Stokes  Theory 

Recently,  progress  has  been  made  in  modelling  the  propagation  of  small 
amplitude  Stokes  waves  due  to  the  simplicity  of  the  dispersion  relation 
incorporating  the  lowost-order  nonlinearity.  Following  Whitham  (1967),  this 
relation  may  be  written  as 

w ^  =  w  1  +  (ka)“D(kh))  (2.1) 

who  re 

2 

=  gktanhkh  (2.2) 

o 

is  the  linear  dispersion  relation,  and  where 


D(kh)  = 


cosh(4kh)  4-8-2  tanh~(khl 
8  sinh^(kh) 


(2.  3) 


(2.1)  may  be  used  to  determine  either  j  or  k  given  a,  the  local  wave 
amplitude,  and  h.  Yue  and  Mel  (1980)  demonstrated  that  (2.3)  may  be 
incorporated  in  the  lowest  order  parabolic  approximation  of  the  Helmholtz 
equation,  yielding 


2ik  A  +  A 
ox  yy 


D  j  A  j  “A  =  0 


(2.4) 


where  A  is  the  complex  amplitude  of  a  wave  described  by 


P  =  Re{ Ae 


i[  k  x  -  u)L  j 

1  n  J  i 


The  resulting  equation  is  in  the  form  of  a  cubic  Schrodinger  equation  for  wave 
evo hit  i  ;M  in  the  x-di recti  in.  This  formulation  has  been  extended  to  the  case 
of  waves  on  currents  in  a  slowly  varying  domain,  yielding  the  equation  (Kirby; 
198),  1984): 


<ci;  *  *  VAy  *  !  h:V^ix  +  (SrlyM 


-  (CCfA  )f  -  a  k  ^  I)  |  A  |  ^  A 


(2.5) 


=  •;  +  kil. 


(2.8) 


">  =  (gktnnhkh) 


(2.7) 


0;  lateral  boundaries 


(2. 28) 


A 

y 

The  computational  model  used  in  the  present  examples  is  obtained  from 
(3.28  -  29)  by  setting  U  =  (U,V)  =  0,  neglecting  any  imposed  ambient 
current.  We  further  adopt  the  lower  order  approximation  of  Radder  (1979)  and 
set  P[  =  0  in  (3.28). 

Data  from  the  laboratory  experiment  of  Berkhoff  et  al  (1982)  are 
available  for  the  labelled  transects  1  through  8  indicated  in  Figure  2.2.  The 
computational  domain  was  discretized  into  square  grids  Ax  =  Ay  =  grid 
spacing),  and  the  grid  scheme  was  established  so  that  grid  rows  coincided  with 
the  measurement  transects.  Grid  size  was  decreased  until  the  point  was 
reached  where  further  reduction  did  not  affect  model  predictions 
significantly.  The  final  numerical  calculations  were  performed  using  a  space 
of  Ax  =  0.25  m.  Results  for  wave  amplitude  normalized  by  the  incident  wave 
amplitude  are  presented  in  Figure  2.3  in  the  form  of  a  contour  plot  for  the 
region  of  the  shoal  and  focus.  The  results  for  the  Hedges  model  are 
qualitatively  similar  to  the  Stokes  waves  results  presented  in  Kirby  and 
Dnlrvnple  (1984),  and  agreement  between  the  Stokes  and  Hedges  results  is 
definitely  better  than  between  either  nonlinear  model  and  linear  theory. 

Plots  if  normalized  amplitude  for  the  labelled  transects  1  -  8  in  Figure 
2.2  i  r  •  ■  given  in  Ki  gur-s  2. 'in  -  h,  respectively.  The  plots  include  results  of 
the  Hedges  model  ,  t  he  Stokes  mode]  ;jnd  the  laboratory  data  of  Berkhoff  et  a  1 
(  1983  ). 
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where  A0  =  0.0232  m  ts  the  amplitude  of  the  incident  wave.  The  wave  period 
T  =  1  sec.  Referring  to  Figure  2.2,  we  establish  slope-oriented  coordinates 

I  » 

(x  ,  y  }  which  are  related  to  the  computational  coordinates  (x,y|  according 
to : 


x  =  (x  -  10. 5)cos  20°  -  (y  -  10)sin  20°  (2.24a) 

t 

=  (x  -  10.  5)sin  20°  4-  (y  -  10)cos  20°  (2.24b) 

If 

The  origin  {x  ,  y  }  =  {0,0}  corresponds  to  the  center  of  the  shoal.  The  slope 
is  described  by: 


0.45  m. 


h  - 


0.45  -  0.02  (5.82  +  x  )m, 


x  <  -  5.82  m 


x  >  -  5.82  m 


(2.25a) 

(2.25b) 


The  boundary  of  the  elliptic  shoal  is  given  by: 

'2  ’2 


(2.26) 


and  the  depth  in  the  shoal  region  is  modified  according  to: 

'  2  '21/2 

h  ‘  h.,t,pe  -  °-5l‘  -  (iff?)  -t1?)]  +  <2’27) 

I  f 

resulting  in  a  depth  h(x  =  0,  y  =  0)  =  0.1332  m. 

The  lateral  boundaries  at  y  =  0,  20  m  are  open,  but  are  far  enough  from 
the  region  of  the  shoal  so  that  we  can  specify  reflective  boundary  conditions 


on  the  lateral  boundaries 
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definitive  laboratory  data  which  shows  the  direct  effect  of  amplitude 
dispersion  in  intermediate  depth.  For  this  reason,  we  include  here  a 
comparison  of  the  refraction-diffraction  calculations  based  on  Hedges  form  of 
the  dispersion  relation,  and  data  taken  from  the  experiments  described  by 
Berkhoff  et  al  (1982).  This  data  set  has  already  been  used  to  check  the 
linear  form  of  the  CREDIZ  model  (Berkhoff  et  al,  1982),  and  a  nonlinear  model 
based  on  the  Stokes  dispersion  relation  (Kirby,  1983;  Kirby  and  Dalrymple, 
1984).  The  reader  is  referred  to  the  previous  literature  for  details  of  the 
experiment.  The  work  of  Kirby  (1983)  and  Kirby  and  Dalrymple  (1984)  compare 
the  Stokes  model  to  results  obtained  using  linear  theory,  and  establishes  that 
the  Stokes  model  provides  a  detailed,  accurate  reproduction  of  the 
experimental  data.  We  therefore  take  the  view  in  the  remainder  of  this 
section  that  the  Stokes  wave  model  represents  the  "correct"  wave  model  for 
waves  in  intermediate  water  depth  (Ur  <  0(1)).  We  will  therefore  be  primarily 
interested  in  how  much  discrepancy  between  numerical  and  experimental  results 
is  introduced  by  using  approximate  dispersion  relations.  We  summarize  the 
important  points  of  the  experimental  arrangement  here.  The  experimental 
topography  consists  of  an  elliptic  shoal  resting  on  a  plane  sloping  bottom 
with  a  slope  of  1:50.  The  plane  slope  rises  from  a  region  of  constant  depth 
h  =  0.45  m,  and  the  entire  slope  is  turned  at  an  angle  of  20°  to  a  straight 
wave  paddle.  Bottom  contours  are  shown  in  Figure  2.2  along  with  the  chosen 
computational  domain,  which  is  indicated  by  the  dashed  line  surrounding  the 
contours.  The  offshore  boundary  of  the  computational  domain  is  chosen  so  that 
water  depth  Ls  constant  along  x  =  0.  The  initial  condition  for  the  wave  then 
correspond;  to  the  uniform  wave  train  generated  at  the  wave  paddle;  we  set: 

A(x  =  O.y)  =  Ao  (2.23) 
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Variation  of  the  Linear  and  shallow  water  approximate  dispersion 

relations  witli  kh  and  wave  steepness  e  =  ka  -  linear; 

*-  '  Hedp.es  (1976);  —  Walker  (1976).  a)  e  =0.1, 

b)  =  0.2,  c)  =  0.3,  d)  e  =  0.4. 


9 

w“/gk  =  (1  +  e/2kh)tanh(kh  +  e/2)  (2.20) 

corresponding  to  Walker's  approximation,  and 

2 

oo  /gk  =  tanh(kh  +  e)  (2.21) 

corresponding  to  Hedges  approximation.  The  variation  is  over  values  of  kh  for 
fixed  values  of  wave  steepness  e  *  ka;  curves  are  given  for  values  of  e  =0.1, 
0.2,  0.3  and  0.4.  The  results  show  that  the  difference  between  Walker's 
approximate  form  and  linear  theory  remains  pronounced  for  large  values  of  kh, 
while  Hedges  approximation  approaches  the  linear  result  move  rapidly.  The 
common  shallow  water  limit  of  the  two  relations  for  small  a/h  is  seen  to  break 
down  for  finite  values  of  a/h,  with  (2.20)  reducing  to 

u)2/  gk  =  kh  +  e  +  e^/4kh  (2.22) 

where  the  last  term  becomes  large  for  any  fixed  e  as  kh  0.  It  seems  that 
Booij's  conclusion  that  Hedges'  approximation  is  preferable  to  Walker's  based 
on  the  deep  water  behavior  extends  as  well  to  the  shallow  water  limit,  due  to 
the  0(e  )  term  which  cannot  be  made  arbitrarily  small  in  shallow  water. 

Hedges  (1976)  approximation  has  been  used  to  form  the  basis  of  nonlinear 
effects  in  at  least  one  operational  refraction-diffraction  model;  the  model 
CREDLZ  developed  by  the  Delft  Hydraulics  Laboratory  in  conjunction  with  the 
Ri jkswaterstaat .  A  recent  calibration  of  this  model  is  described  in 
Dingemanns  (1983)  and  Dingemanns  et  al  (1984).  To  date,  no  adequate 
comparison  has  been  made  between  the  approximate  nonlinear  model  and  any 


9 


Alternat ive Ly ,  Hedges  (1976)  proposed  a  similar  relation  of  the  form 


2 

co 


gk  tanh(k[h  +  a  ]  ) 


(2.16) 


Both  (2.15)  and  (2.16)  have  similar  properties  in  the  limits  of  deep  and 
shallow  water.  In  shallow  water  and  for  a/h  small,  both  formulas  lead  to  a 
dispersion  relation 


2 

u> 


gk2(h 


+  a) 


or,  equivalently,  a  phase  speed 


(2.17) 


C  =  [g(h+a))l/2  (2.18) 

This  phase  speed  represents  the  speed  of  propagation  of  a  solitary  wave  of 
height  H  =  2a.  The  shallow  limit  of  either  approximate  form  is  thus 
physically  reasonable. 

In  the  deep  water  limit,  both  (2.15)  and  (2.16)  approach  the  linear 
dispersion  relation,  as  the  ratio  a/h  approaches  zero  due  to  increasing  water 
depth.  Booij  (1981)  has  mentioned,  without  demonstration,  that  the  approach 
to  the  linear  form  is  quite  slow  if  the  relation  (2.15)  of  Walker  is  used. 
Figure  2.1  shows  the  variation  of  the  right-hand  side  of  the  relations 

2 

— r—  =  tanhkh  (2.19) 


cor  respond ing  to  linear  theory, 
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Recently,  Liu  et  al  (1985)  have  shown  that  refraction-diffraction 


calculations  may  be  performed  for  regular  cnoidal  waves  using  a  spectral 
approach  based  on  the  Boussinesq  equation.  This  method  is  reasonably  simple 
to  develop,  and  computations  are  straight  forward  and  fairly  economical  due  to 
the  simplicity  of  the  wave-wave  interaction  coefficients  in  shallow  water. 
However,  the  calculations  are  highly  dependent  on  the  relative  phase  speeds 
and  amplitudes  of  a  number  of  waves  in  a  discrete  spectrum,  and  the  method  can 
not  be  applied  to  the  present  monochromatic  model. 

Alternat ivelv,  several  authors  have  proposed  simple  modifications  to  the 
linear  dispersion  relation  which  are  designed  to  mimic  the  effect  of  amplitude 
dispersion  in  shallow  water.  Based  on  a  large  number  of  laboratory 
observations  of  broken  and  unbroken  waves  propagating  over  a  focusing 
topography,  Walker  (1976)  proposed  that  the  nonlinear  effects  in  shallow  water 
could  be  modelled  in  a  refraction  scheme  by  modifying  the  predicted  linear 
phase  speed  according  to 


Cn!l  =  V1  +  a/h)  (2.13) 

where  "a"  is  the  wave  amplitude  and 

CZ  =  k=  (gtanhkh/k)1/-  *  (2.14) 

Bool  i  (1981)  showed  that  a  dispersion  relation  equivalent  to  (2.13)  may  be 
written  as 


J  =  1  +  tip  tanh(k[h  +  §"i) 


(2.15) 


The  limitations  of  the  Stokes  theory  leads  to  two  modelling  problems: 

1)  An  alternate  theory  (cnoidal  waves)  must  be  used  to  describe  waves  in  the 
range  where 

iAl/h  „  0(1); 

(kh; 

2)  A  way  must  be  found  to  describe  the  smooth  transition  of  waves  from  the 
Stokes  regime  into  the  cnoidal  regime. 

The  first  problem  results  from  the  limits  of  validity  of  the  Stokes 
theory.  The  second  problem  arises  due  to  the  fact  that  wave  properties  do  not 
vary  smoothly  across  the  division  between  Stokes  and  cnoidal  waves. 

Using  the  results  of  cnoidal  wave  theory,  the  dispersion  relation  for 
periodic  waves  may  be  written  as  (Flick,  1978) 

C  =  (gh) l/*2  (  1  +  f  j(m))  (2.11) 

whe  re 

fl(m)  =  m  ^ 2  "THO-  (2*l2) 

and  where  E(m)  and  K(m)  are  the  Jacobean  elliptic  integrals. 

This  relation  may  be  used  to  calculate  local  values  of  k  in  a  modelling 

scheme.  This  analytic  approach  has  been  used  in  several  refraction  schemes 

(Chu,  1975;  Skovgaard  and  Petersen,  1977;  Headland  and  Chu,  1984)  as  a  means 

for  Including  the  effect  of  shallow  water  amplitude  dispersion.  The  full 

formulation  requires  a  large  amount  of  iterative  computation  due  to  the 

interrelationship  of  the  wave  height  H  and  elliptic  parameter  m. 
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Here,  w  is  the  fixed  wave  frequency  with  respect  to  the  stationary  domain,  and 
the  wave  is  assumed  to  be  propagating  in  the  x-direction. 

Equations  of  the  form  (2.5)  (or  the  large  angle,  large  current  extension 
described  in  Chapter  3)  are  applicable  in  intermediate  depth  ( kh  ~  0(1))  and 
are  correct  in  the  deep  water  asymptote  (kh  +  “)  under  the  condition  of 
k | A |  <<  l.  In  shallow  water,  the  regular  expansion  in  k | A |  breaks  down  and  is 
replaced  instead  by  an  expansion  to  lowest  order  in  kh  «  1,  |A|/h  «  1,  with 
the  Ursell  number 


U 

r 


1  A1  /h 

(kh)2 


0(1) 


(2.8) 


representing  the  region  where  cnoidal  and  solitary  waves  are  described.  The 
breakdown  in  the  Stokes  theory  may  be  seen  by  looking  at  the  shallow  limit  of 
(2.3); 


D(kh) 


kh  +  0 


(2.9) 


The  singularity  in  D  as  kh  +  0  is  quite  severe  and  overwhelms  the  computations 
if  regions  of  the  computational  domain  become  too  shallow.  The  dispersion 
relation  may  be  arranged  to  the  form 


2 

u) 


(2.10) 


indicating  that  | A | /h  must  remain  extremely  small  as  kh  *  0,  in  order  for  the 
Stokes  theory  to  remain  valid.  This  requirement  is  clearly  not  applicable  in 
the  vicinity  of  surfzones,  where  | A | /h  is  typically  of  0(1). 
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Figure  2.4. 


Comparison  oE  model  using  Hedges  (1976)  dispersion 
relation  to  the  Stokes  model  (Kirby  and  Dalrymple, 
1984)  and  laboratory  data  (Berkhoff  et  al . ,  1982). 

-  Stokes  model;  -  Hedges  model;  o  laboratory 

data,  a-h)  sections  1-8,  respectively. 


Matching  Stokes  and  Shallow  Water  Dispersion 


The  inadequacy  of  an  approximate  shallow  water  dispersion  to  model 
nonlinear  effects  in  intermediate  water  depth,  coupled  with  the  invalidity  of 
the  Stokes  model  in  shallow  water,  leads  to  the  need  for  a  matched  dispersion 
relation  which  predicts  the  phase  speeds  of  waves  smoothly  from  deep  to 
shallow  water.  This  problem  is  central  to  the  prediction  of  the  properties  of 
shoaling  waves. 

Goring  (1978)  has  shown  that  a  matching  exists  between  Stokes  and  cnoidal 
waves  in  the  limit  of  small  amplitude,  by  showing  that  a  regular  perturbation 
solution  of  the  Korteweg-deVr ies  equation  is  equivalent  to  the  shallow  water 
limit  of  the  Stokes  solution  for  general  water  depth.  However,  for  larger 
amplitudes  the  series  solution  does  not  converge,  and  the  solution  in  terms  of 
powers  of  c  =  ka  becomes  inappropriate.  In  effect,  as  water  depth  decreases, 
the  Stokes  solution  may  become  invalid  before  the  region  of  validity  of  the 
cnoidal  theory  (small  kh)  is  reached.  Flick  (1978)  suggested  that  it  may  be 
possible  to  construct  a  dispersion  relation  by  means  of  matched  asymptotic 
expansions;  he  suggests  a  relation  of  the  form 

C  =  r=C  +  C  -  C.  .  (2.29) 

k  s  c  £in 

where 

2  1/2  1/2 

Cg  =  ( 1  +  e  D)  (gtanhkh/k)  (2.30) 

is  the  Stokes  phase  speed, 


C  =  (gh)1/2(l  +  ”  f.(m)) 


f  L  (m)  =  (2  -  3E/m)/K(m)  -  m)/m  (2-31) 

i 

i 

is  the  cnoidal  phase  speed,  and 

CMn  =  <Sh)l/2(  1  "  ^<kh>2  +  (2.32) 

(kh)2 

is  both  the  shallow  water  limit  of  (2.30)  and  the  small  amplitude 

(m  *  0)  limit  of  (2.31).  The  fact  that  these  limits  coincide  lends  further 

creedence  to  the  conclusions  of  Goring  (1978).  However,  in  cases  where  the 

small-amplitude  limit  is  inappropriate,  (2.32)  may  not  be  used  directly  as  a 

reasonable  approximation.  This  would  seemingly  include  the  majority  of 

2 

cnoidal  wave  cases,  where  the  Ursell  parameter  (a/h)/(kh)  is  0(1). 

For  the  purpose  of  inclusion  in  the  monochromatic  wave  model,  we  propose 
that  a  dispersion  relation  of  the  form 

w2  =  gk(  1  +  f  |  (kh)t  ~ d]  tanh(kh  +  f2(kh)e)  (2.33) 

may  be  constructed  In  order  to  model  nonlinear  effects  over  a  broad  range  of 
depths.  Comparing  (2.33)  to  the  previous  forms  of  the  dispersion  relation,  wi 
see  that  the  Stokes  wave  form  Is  recovered  by  the  choice 


f  ^ (kh)  =  1. 


all  kh 


(2.34) 


f  2 ( kh )  =  0. 


while  the  Hedges  model  is  recovered  by  choosing 


0 


fj(kh)  = 


all  kh 


(2.35) 


f  2 ( kh )  =  1. 


The  composite  model  may  be  constructed  by  choosing  forms  for  the  arbitrary 
functions  f^  and  f2*  In  particular,  we  require  that  f^(kh)  ■+■  1  as  kh  +  00 
while  f2(kh)  ■*  0  as  kh  >  00 ,  in  order  to  recover  the  Stokes  wave  limit.  In 
shallow  water,  we  require  f,  ♦  1  as  kh  *  0,  while  f^(kh)  must  be  of  0(kh"*)  or 
smaller  in  order  to  overcome  the  singularity  in  D,  which  is  0(kh  ).  Based  on 
these  requirements,  we  have  chosen  f.  and  t‘2  according  to 


f^(kh)  =  tanh  kh 


(2.36) 


f„(kh)  =  [ kh/sinh(kh)] 


(2.37) 


The  dispersion  relation  resulting  from  these  choices  is  illustrated  in 
Figure  2.5,  where  the  right-hand  sides  of  the  following  relations  are  plotted. 


Linear 

2 

; —  =  tanhkh 


(2.38) 


Stokes 
(  2 

=  (1  +  e  D)tanhkh 


(2.39) 


Hedges 

2 

^7  =  tanh(kh  +  e ) 


(2.40) 
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Composite 

2  2 

=  (1  +  f.e  D)  tanh(kh  +  f„e) 

gk  1  2 


(2.41) 


As  in  Figure  2.1,  results  are  plotted  for  a  range  of  kh  values,  and  separate 
sets  of  curves  are  given  for  choices  of  e  =  0.1,  0.2,  0.3  and  0.4.  An 
inspection  of  each  set  of  curves  shows  that  the  composite  form  (2.41)  matches 
smoothly  to  both  the  Stokes  form  (2.39)  in  deeper  water  and  the  Hedges  form 
(2.40)  in  shallow  water.  The  strong  singularity  in  the  Stokes  form  is  also 
apparent.  The  curve  corresponding  to  the  composite  form  typically  lies 
between  the  Stokes  curve  and  the  Hedges  curve,  which  may  or  may  not  intersect 
each  other. 

We  remark  that  the  choice  of  the  exponent  4  in  (2.37)  led  to  the  best 
positioning  of  the  composite  curve  (2.41)  with  respect  to  both  the  Stokes  and 
Hedges  relations,  with  both  higher  or  lower  exponents  leading  to  fairly  large 
deviations  over  some  portion  of  the  range  of  kh.  The  choice  of  the  exponent  5 
in  (2.36)  is  of  course  dictated  by  the  requirement  imposed  by  the  singularity 
in  D. 

The  composite  model  was  tested  by  recomputing  the  example  of  Berkhoff 
et  al  (1982)  using  the  modified  dispersion  relation.  The  dispersion  relation 
is  incorporated  directly  in  the  nonlinear  Schrodinger  equation  for  complex 
amplitude  A  following  the  method  of  Kirby  and  Dalrymple  (1984).  The  nonlinear 
modification  to  the  time-dependent  form  of  the  linear  mild-slope  equation  for 
velocity  potential  >>  at  the  free  surface  is  given  by 


P  -  7  •  (CC  7  $)  +  (w2  -  k2CC  )p  +  Ft  =  0  (2.42) 

tt  n  g  n  g 


whore 


„  2  2 
F  =w  -  w 

o 


(2.43) 


2 

and  is  the  linear  wave  frequency  squared; 


-  gktanhkh 


(2.44) 


2  ~ 

Representing  w  by  (2.33)  and  letting  $  be  represented  as 


.  l(k  x  -  u  t) 

=  Re{ -  A(x,v)e  ° 


(2.45) 


where  k  is  a  reference  wave  number,  reduces  (2.42)  to  the  form 


C  A  +  1C  (k  -  k)A  +  -kc  )  A  -  V,  •  (CC  V  A) 
2  x  2  0  2  r  2o)  h  2  h 

x  o 


Icj 


+  2  Ltl  +  fik  lAl  D) 


„  0  tanh(kh  +•  f  ,.k  |  A|  ) 


tanh  kh 


-  l]  A  =0  (2.46) 


This  form  of  the  nonlinear  model  is  used  to  compute  results  corresponding  to 
the  experimental  data.  The  computational  grid  and  domain  is  as  described  in 
section  2.2.  A  parabolic  form  of  (2.46)  is  obtained  by  using  the  scaling 
assumption 


|Aj  «  0(  ik  |  A  |  )  ,  (2.47) 

leading  to  the  neglect  of  the  term  (CC  A  )  in  (2.46).  Again,  the 

O  XX 

computational  scheme  is  described  in  Chapter  3;  here  we  use  the  choice  Pj  =  0 
A  contour  plot  of  normalized  amplitude  [ A | / A0  is  given  in  Figure  2.6. 

The  contours  resul  ing  from  calculations  using  the  composite  model  are  closer 
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guru  2. 


Ampl i tude  contours  relative  to  incident  wave  amplitude 
present  composite  model. 
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in  shape  to  the  results  of  the  Stokes  model  than  are  the  contours  of  the 
Hedges  model,  as  would  be  expected.  Plots  of  normalized  amplitude  along  the 
labelled  transects  1-8  are  presented  in  Figure  2.7,  where  the  composite 
model  results  are  compared  to  Stokes  model  results  and  laboratory  data.  The 
results  of  the  composite  model  are  generally  closer  to  the  Stokes  model  and  to 
the  data  than  are  the  results  of  the  Hedges  model,  although  the  discrepancies 
in  either  case  are  not  large,  with  all  models  generally  showing  better 
agreement  with  data  than  the  linear  model. 

2.4  Conclusions 

The  present  chapter  has  presented  a  method  for  extending  the  effects  of 
nonlinear  dispersion  in  a  monochromatic  Stokes  wave  model  into  water  depths 
which  are  too  shallow  for  Stokes  theory  to  retain  its  validity.  The  proposed 
model  provides  a  smooth  patch  between  Stokes  theory  and  an  empirical  shallow 
water  relation  due  to  Hedges,  with  the  two  separate  forms  being  obtained  in 
the  limit  of  deep  and  shallow  water,  respectively. 

Although  example  calculation  show  that  the  differences  between  the  Stokes 
model  and  Hedges  model  are  not  large  In  intermediate  water  depths,  we  feel 
that  the  present  results  are  preferable  to  the  use  of  the  Hedges  relation 
alone  for  several  reasons.  First,  the  fact  that  the  Hedges  model  approaches  a 
linear  model  In  deep  water  implies  that  discrepancies  would  inevitably  exist 
in  calculations  which  start  or  extend  into  fairly  deep  water.  Existing 
laboratory  data  covers  too  shallow  a  depth  range  (and  too  low  )  range  of  wave 
heights)  to  show  these  discrepancies  fully.  The  present  method  provides  a 
moans  for  patching  the  empirics 1  sttallow  water  model  into  a  more  reasonable 
deep  wat er  I i mi t . 


Secondly,  although  the  difference  between  results  of  the  Hedges  and 
Stokes  models  Is  not  severe,  the  errors  In  phase  and  group  velocity  would  be 
expected  to  accumulate  over  long  distances  for  waves  propagating  in  shallow 
water.  The  proposed  composite  model  reduces  the  small  local  error  in  the 
Hedges  model  by  mixing  the  Hedges  and  Stokes  effects  together,  and  would  be 
expected  to  produce  more  accurate  results  for  waves  propagating  over  large 
distances. 
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The  remaining  term  C4  is  given  by 


w  .  io 
2  A  +  2 


a  1  +  f1(k|A|)2D) 


tanh(kh  +  f 7k| A| ) 
tanh  kh 


which  is  written  out  in  finite  difference  form  as 
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The  computations  proceed  by  updating  values  of  A  from  the  known  "i"  row  to  the 
unknown  ''i+l"  row.  Iterations  for  the  nonlinear  terms  are  performed  using  a 
repeated  implicit  calculation  as  in  Kirby  (1983). 

Using  the  notation  (3.30),  (3.28)  is  written  in  finite  difference  form  as 
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where  Cl,  C2,  C3,  CPI,  CP2  and  CP3  are  known  coefficients  given  by 
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(3.28) 


This  equation  extends  the  published  results  in  Kirby  (1984)  and  Kirby  and 
Dalrymple  (1983b)  to  the  case  of  waves  propagating  at  large  angles  to  the 
x-direction  in  the  presence  of  a  current.  Nonlinearity  here  is  consistent 
with  the  lowest  order  Stokes  formulation  of  Kirby  and  Dalrymple  (1983a).  The 
nonlinear  term  may  be  directly  modified  according  to  the  results  in  Chapter 
2.  The  resulting  equation  is  of  the  form 

tanh(kh  +  f  k | A | ) 

(C  +  U)A  +  ~ {  (1  +  f .  (k|A|  )  D  - ; - ~r -  -  l}  A  +  ...  =  0  (3.29) 

g  x  2  I  *  1  tanh  kh 

where  fj,  D  are  functions  of  kh  given  by  (2.36),  (2.  37),  and  (2.3  ) 

respectively. 


3. 2  Finite  Difference  Scheme 

The  finite  difference  scheme  for  (3.28)  follows  directly  using  the  Crank- 
Nlcolson  method  for  performing  an  implicit  update  for  each  row  in  x.  The 
notation  of  Kirby  (1983),  section  6.2  is  retained.  We  denote  x-positions  by 
”t"  superscripts  and  '/-positions  by  -’j”  subscripts.  For  example: 
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where  a  and  a^  are  a^  and  a^  written  in  terms  of  A.  The  first  bracketed 
term  is  given  by 

{ Socj  '  +  SCR^e-1*  +  a 3 ' } 

=  3{2it)U(^)x+  i(o>U)x(^)  -  2WUk(e  -  l)(f) 

+  2iaVC^)y  +  i(°v)y(^  "  i(k<e  “ 

-  2ik(e  -  l)UV(£)y  -  (uv(£)x)y  -  (Uv(£)y)x 
+  [<P  -  v2)  $y]y} 

*t  <  <“'%(£)  *  <“”>„$  +  (“V>y^)x+  ik(E 

+  3<»U>x(j-)x  +  3ik(e  -  1)(»U)X  J-  1 

So  far,  nothing  has  been  dropped  besides  some  second  derivatives  of 
currents.  In  order  to  simplify  the  equation  further,  we  drop  the  following 
terms  as  being  small: 

1)  Terms  in  (1  -  e)  multiplying  products  of  current  components  or 
derivatives . 

2)  Terms  with  8  multiplying  derivatives  of  the  ambient  current. 

We  are  then  left  with  the  form  of  the  equation  to  be  used  in  the  numerical 
scheme. 


-  i>(«v)y$ 

(3.27) 
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P2  -  P^k/k) 


p2  -  pi  ~  k/k) 


=  Pt(l  -£) 


(3.25) 


Denote  £  =  k/k.  It  is  usually  assumed  that  e  =  1;  i.e.,  the  wave  number 

varies  only  slightly  along  a  grid  row.  Making  this  assumption,  we  may  take 
1  -  e  «  0(1),  introducing  a  small  parameter.  The  full  equation  for  A  defined 
according  to 
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The  ri.RhL-hand  side  of  (3.17)  becomes 
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The  R.H.S.  of  (3.17)  reduces  to 
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We  now  need  to  simplify  Rj  and  R?.  Shifting  to  the  reference  phase  ^  defined 
in  (3.6),  we  continue  to  expand  all  terms  in  order  to  eliminate  the  phase 
function  v  and  >>.  The  use  of  i3  injects  a  k  in  the  equation  in  place  of  k  in 


some  terms.  Note  that 
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will  have  some  contributions  like 


ri  i 
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These  will  balance  the  {  }  terms  above  in  (3.16)  and  give  the  lowest  order 
diffraction  and  y-direction  transport  terms.  Terms  without  y  derivatives  may 

f  ~ 

he  manipulated  simply.  In  terms  of  the  complex  amplitude  A  ,  M<J>  becomes 
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where  the  phase  is  successfully  eliminated  from  dp  Differentiating  (3.19) 
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(3.13)  is  then  given  by 
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The  term  {  }  will  be  cancelled  by  leading  order  contributions  from 


Substituting  for  <j)  using  (3.11)  gives 
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The  right-hand  side  remains  to  be  simplified.  Note  that 
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Finally,  a  reference  phase  function  is  defined  according  to 
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is  the  average  of  the  linear  wave  number  over  a  y-coordinate  grid  row.  We  now 

proceed  by  reducing  (3.1)  to  a  usable  form. 

2 

Denote  k(p  -  U  )  =  y.  Then,  (3.1)  may  be  written  as 
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gradient  operator,  and  w  represents  a  damping  coefficient  due  to  the  presence 
of  boundary  layers  or  porous  bottom  damping  (Dalrymple  et  al,  1984)  or  wave 
breaking  (Kirby,  1983).  Kirby  (1983)  developed  (3.1)  with  U  =  0,  and  Pj  =  1/4 
in  its  full  form  and  gave  the  numerical  scheme  for  the  large  angle 
approximation  without  currents.  This  chapter  provides  the  same  information 
for  the  case  of  waves  interacting  with  a  pre-speci f ied  ambient  current. 

Section  3.1  gives  the  derivation  of  an  evolution  equation  for  complex 
amplitude  A  based  on  (3.1)  with  Pj  =  1/4.  Section  3.2  provides  the  finite 
difference  form  of  the  evolution  equation  using  the  Crank-Nicolson 
approximation. 

3.1.  The  Parabolic  Equation  for  Waves  on  Currents 

As  a  first  step,  the  parabolic  equation  for  4>  is  developed.  Several 
simplifying  assumptions  will  be  employed.  First,  nonlinear  terms  and 
dissipative  terms  will  be  considered  to  be  grouped  with  the  highest  order 
contributions  in  the  equation,  and  only  the  largest  contributions  resulting 
from  manipulations  to  these  terms  will  be  retained.  For  example: 

(o2k2D|A|2  $)  *  a2k2D|A|2  »  ia2k3DiA|2  £  (3.5) 

dx  '  ax 

These  terms  will  therefore  not  retain  any  derivatives  related  Lo  the  slow 
spatial  scale  of  variation. 

Booij  (1981)  drops  all  terms  containing  products  of  ambient  velocity 
components.  Here,  these  terms  will  he  retained  up  to  the  level  of  second 
derivatives;  i.e.,  terms  such  as 
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Chapter  3.  Large  Angle  Formulation  for  the  Wave-Current  Model 


Kirby  (1983)  provided  a  higher-order,  large  angle  formulation  of  the 
parabolic  approximation  for  waves  on  currents,  following  the  original  work  of 
Booij  (1981).  The  governing  equation  for  the  value  of  the  velocity  potential 
$  at  the  mean  water  level  is  given  by  Kirby  (1983)  as 


where 
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The  choice  Pj  =  0  reduces  the  model  to  the  lowest  order  approximation  of 
Radder  (1979)  or  Kirby  (1983,1984).  The  choice  Pj  =  1/4  gives  the  higher 
order  approximation  proposed  by  Rooi j  (1981).  Finally,  U  =  (U,V)  represents 
the  ambient  current  In  x  am.  y  coordinates,  Vh  represents  the  horizontal 


r< 


~  i+ 1  i+  1 

Here,  the  unknown  A.  is  given  by  the  last  updated  value  of  A.  in  the 

J  J 

iteration  procedure.  Chapter  4  describes  how  wave  breaking  is  incorporated  in 
the  specification  of  the  dissipation  coefficient  w. 


t 

i 
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Chapter  4.  Modelling  Waves  In  Surfzones  and  Around  Islands 


4.  1.  Introduction 

In  recent  years,  the  parabolic  equation  method  for  wave  propagation  has 
seen  a  rapid  development  in  the  context  of  predicting  surface  water  waves  over 
areas  of  variable  bathymetry,  including  both  the  effects  of  refraction  and 
diffraction.  The  original  development  of  the  model  was  based  either  on  a 
splitting  of  the  elliptic  mild  slope  equation  of  Berkhoff  (1972)  into  coupled 
equations  for  forward-  and  back-scattered  wave  motion,  as  in  Radder  (1979),  or 
on  direct  perturbation  expansion  of  the  governing  equations  using  the  WKB 
formalism,  as  in  Liu  and  Mei  (1976).  Recent  extensions  of  the  capabilities  of 
the  parabolic  method  include  the  modelling  of  wave-current  interaction  (Booij 
(1931),  Kirby  (1984)),  iterative  calculation  of  the  reflected  wave  field  (Liu 
and  Tsay  (1983)),  and  the  inclusion  of  lowest  order  nonlinear  effects  in  the 
Stokes  wave  formulation  (Kirby  and  Dalrymple  (1983,  1984),  Liu  and  Tsay 
(1984)). 

An  advantage  of  the  parabolic  method  over  solution  techniques  for 
elliptic  and  hyperbolic  equations  is  that  no  downwave  boundary  condition  is 
needed  for  the  solution  of  the  initial  boundary  value  problem.  However,  in 
applications  of  wave  models  to  coastal  areas,  accurate  modelling  of  the 
behavior  of  waves  in  the  vicinity  of  a  physical  downwave  boundary  consisting 
of  an  actual  coastline  or  an  offshore  island  is  of  primary  importance  to  the 
prediction  of  known  physical  effects  such  as  wave-induced  runup,  longshore 
currents  and  scour. 

Wave  breaking  in  the  surf  zone  is  a  complex,  highly  nonlinear 
phenomenon.  It  is  obvious  that  the  parabolic  equation  method,  which  is 
limited  to  the  representation  of  linear  or  weakly  nonlinear  wave  fields,  is 
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basically  Incapable  at  present  of  representing  the  underlying  physics  of  the 
breaking  process.  However,  some  progress  can  be  made  by  shifting  our  view  of 
the  model  from  its  physical  basis  to  its  use  as  a  predictive  tool. 

The  forces  leading  to  the  generation  and  maintenance  of  setup  and  wave- 
induced  currents  depend  on  a  physical  balance  between  gradients  of  excess 
momentum  fluxes,  pressure  forces  due  to  changes  in  mean  surface  elevation,  and 
bottom  shear  stresses.  The  role  of  a  wave  model  in  determining  the  balance 
consists  of  predicting  the  local  wave  energy  density  and  direction  of 
propagation  of  the  wave  field.  Thus,  as  a  lowest  approximation  of  the  overall 
physics,  it  suffices  that  the  wave  model  be  able  to  predict  the  local  wave 
amplitude  in  the  breaking  zone  with  some  degree  of  reliability. 

The  simplest  model  of  wave  decay  in  the  surf  zone,  the  spilling  breaker 
model,  is  based  on  the  assumption  that  the  ratio  of  wave  height  to  local  water 
depth  has  the  same  value  everywhere  in  the  surf  zone  as  at  the  breaker  line. 
This  assumption  has  been  used  extensively  in  the  literature,  from  predictions 
of  setup  (Longuet-Higgins  and  Stewart  (1963))  and  longshore  currents  (Longuet- 
Higgins  (1970))  up  to  the  latest  applications  of  numerical  refraction  schemes 
to  the  study  of  wave-induced  circulation  over  arbitrary  bottoms  (e.g. , 

Ebersole  and  Dalrymple  (1980)).  However,  it  has  long  been  known  that  breaking 
waves,  especially  of  the  plunging  type,  do  not  follow  so  simple  a  rule. 

E  .tensive  model  tests  of  normally  incident  wave  trains  breaking  on  laboratory 
beaches  have  shown  that  the  pattern  of  wave  height  decay  across  the  surf  zone 
is  strongly  a  function  of  the  beach  slope.  Representative  measurements  of 
Horikawa  and  Kuo  (  1966)  are  shown  for  example  in  Figure  A.  I  in  comparison  to 
their  dissipation  model. 

The  purpose  of  the  present  study  is  to  relate  an  empirical  model  of  surf 
zone  wave  energy  decay  to  the  dissipation  coefficient  w  of  the  dissipative 
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wave  model  of  Dalrymple,  Kirby  and  Hwang  (1984),  and  to  detail  the  application 
of  the  model  to  the  prediction  of  wave  height  in  the  surf  zone.  Here,  the 
model  of  Dally,  Dean  and  Dalrymple  (1984)  is  used,  although  any  of  the  related 
models  for  dissipation  in  bores  could  be  applied  just  as  well  (Horikawa  and 
Kuo,  1966;  Divoky  et  al ,  1970;  Battjes  and  Janssen,  1978,  for  example). 

4.2.  The  Energy  Decay  Model 

Dally,  Dean  and  Dalrymple  (1984)  have  proposed  that  the  decay  of  energy 
flux  with  distance  in  the  surf  zone  is  proportional  to  the  excess  of  energy 
flux  over  a  stable  value,  given  for  waves  propagating  shoreward  in  the  x 
direction  by  the  relation 

(EC  )  =  -  ~  (EC  -  (EC)J  (4.1) 

3x  g  h  *•  g  g  s' 

where  h  is  the  local  water  depth  and  K  is  a  constant  to  be  determined,  which 

is  related  to  the  rate  of  energy  decay.  The  quantity  (ECg)s  is  the  "stable" 

2 

energy  flux  for  a  broken  wave  in  water  of  depth  h.  Here  E  is  1/8  pgH  ,  H  is 

2kh 

the  local  wave  height  and  Cijr  =  (l  +  stnh  2kh^  ^  w*iere  p  Is  t'le  fluid  density, 
g  is  the  acceleration  of  gravity  and  k  and  h  are  related  by  the  dispersion 
relationship,  oo  =  gk  tanh  kh.  Here  w  =  2ir/T,  where  T  is  the  wave  period. 
Dally  et  al .  show  that  this  model  of  wave  energy  decay  is  analogous  to  the 
energy  loss  in  a  hydraulic  jump.  The  stable  energy  flux  may  be  related  to  the 
height  obtained  asymptotically  by  a  wave  propagating  over  a  flat  bottom  or  a 
plane  slope.  Measurements  by  Horikawa  and  Kuo  (1966)  indicate  that  a  value  of 
lig  -  0.4  h  is  approached  for  waves  breaking  on  a  plane  slope.  In  the 
following  derivation,  we  denote  (H/h)  =  y,  where  Y  is  an  empirical  constant. 


Equation  (4.1)  may  be  related  to  a  wave  energy  equation  for  dissipative 
motion  after  assuming  a  time-steady  wave  field.  For  normally  incident  waves, 
the  energy  equation  becomes 


r—  (EC  )  =  -  WE 
3x  g 


(4.2) 


Noting  that  C 

gs 


Cff,  W  may  be  written  as 

O 


w 


KC  E  K  C 


(4.3) 


where  H  =  2|a|,  and  A  Is  a  generally  complex  measure  of  wave  amplitude  and 
phase. 

For  a  plane  beach  with  slope  s,  and  neglecting  the  effect  of  setup  (which 
is  not  calculated  by  the  wave  model  directly),  a  simple  analytic  solution  of 
(4.1)  was  obtained  by  Dally  et  al.  Their  comparison  to  laboratory  data  of 
breaking  waves  shows  that  this  model  very  successfully  represents  wave  height 
decay  across  the  surf  zone.  This  analytical  model  will  therefore  be  used  for 
comparison  with  the  numerical  model.  Defining  a  =  K/s,  the  wave  height  in 
dimensionless  form  is  related  to  the  local  water  depth  by: 


(ir)  -(£-)  [ci^Xfrd 

b  b  b 


o-5/2 


(4.4) 


where 
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breaking) 
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tj-)  =  (h~)  t  1  ’I®  5  a  =  y  •  (4.6) 

b  b  b 

Based  on  comparisons  Co  Che  laboratory  data  of  Horikawa  and  Kuo  (1966), 
Dally  et  a  1.  chose  the  value  K  =  0.15.  The  special  case  a  =  5/2  then 
corresponds  to  a  beach  slope  s  =  0.06.  Results  for  a  range  of  a  values  of 
l  <  a  <  10,  corresponding  Co  Che  range  of  beach  slopes  0. 15  >  s  >  0.015,  are 
given  In  Figure  4.2  for  y  =  0.4,  k  =  0.78.  The  lines  labelled  1  and  2 
correspond  to  the  constant  decay  H  =  xh  =  0. 78  h  and  to  the  stable  wave  height 
Hs  =  yh  =  0.4  h,  respectively.  The  effect  of  beach  slope,  s,  is  clearly 
apparent.  For  mild  slopes,  a  is  large  and  the  wave  heights  across  the  surf 
zone  are  much  less  than  given  by  the  spilling  breaker  model,  while,  for  steep 
beaches,  a  is  small  and  the  wave  heights  exceed  the  spilling  breaker  model 
(H  =  <h)  heights.  The  paper  by  Dally  et  al .  contains  an  extensive  comparison 
between  these  results  and  the  data  of  Horikawa  and  Kuo  (1966). 

4.3.  Application  of  the  Model  in  the  Parabolic  Equation  Method 

The  results  (4.4-6)  provide  a  check  for  determining  the  accuracy  of 
iterative  schemes  using  the  wave  damping  term  (4.3).  Noting  that  Hg  =  yh  and 
H  =  2 | A | ,  (4.3)  may  be  written  as 

KC  2,2 

M  — [  1  -  ^  •  (4.7) 

‘l*l! 

Assuming  that  the  reflected  wave  (in  the  minus  x  direction)  is  small,  we  use 
the  parabolic  equation  given  by  (3.28)  to  obtain 

Ay  "  i(k-k  )A  +  — ^ —  C  A  -  ~ —  (CC  A  )  +  wA  =  0  (4.8) 

x  o  2C  gx  2C  g  y  y 

8  8 


where  we  have  neglected  currents  and  nonlinear  effects. 


Subscripts  x  and  y  denote  differentiation.  Here,  A(x,y)  is  the  complex 
amplitude  of  a  steady  wave  train 

i(k  x  -  ut) 

n(x,y,t)  =  Re{  A(x,y)  e  °  }  ,  (4.9) 


kQ  is  a  reference  value  of  the  wavenumber  based  on  a  constant  water  depth 
which  is  characteristic  of  the  domain,  and 


w 


JL  =  —  f  i 

2C  2h  1 


g 


2  2 
y  h 


4  A 
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Here  we  choose  a  real  value  of  w  in  contrast  to  the  results  for  boundary  layer 
damping  (Liu,  1984),  since  the  wave  breaking  process  would  not  be  expected  to 
distort  the  wavelength  in  the  same  manner  as  a  small  dissipation  due  to 
viscous  effects.  Equation  (4.8)  is  written  in  finite  difference  form 
according  to  the  Crank-Nicolson  method: 
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whore  v-derlvative  terms  are  given  by 
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The  intermediate  value  A.  is  provided  by  information  on  the  starting  grid 
row  according  to 
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7  or  already  breaking 
h  2 


or  not  yet  broken 

h  2 


(4.15) 


~i+ 1 

If  waves  are  unbroken  in  all  j  grid  blocks,  w^  is  set  equal  to  zero  and  the 
scheme  proceeds  based  on  equation  (4. 14)  with  no  second  iteration. 
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4.4.  Nomallv  Incident  Waves  on  a  Plane  Slope 


Several  cases  were  run  for  waves  starting  in  a  depth  of  2  m  and 
propagating  directly  towards  shore  over  a  plane  slope.  The  program  checks  at 
each  step  that  the  wave  height  has  not  exceeded  the  breaking  criterion.  When 
H  becomes  greater  than  <h,  the  program  begins  calculating  values  of  the 
damping  coefficient  (4.15).  Breaking  continues  until  w  falls  to  a  value  of 
zero,  which  does  not  occur  on  the  plane  beaches  studied  here  but  would  be 
expected  to  occur  readily  for  waves  propagating  over  uneven  topography.  For 
each  case  examined,  the  wave  height  was  assumed  to  be  1.0  m  at  a  depth  of  2  m, 
and  wave  period  was  assumed  to  be  5  seconds.  Values  of  a  =  1,  3,  and  10  were 
tested  using  various  computational  grid  spaces  Ax.  Results  for  a  =  1  and 
Ax  =  0.2  m  and  1.0  m  are  shown  in  Figure  4.3,  a  =  3  and  Ax  =  1.0,  2.0,  and 
5.0  m  in  Figure  4.4  and  a  =  10  and  Ax  =  2.0,  5.0  and  10.0  m  in  Figure  4.5. 

The  exact  solution  (4.4)  is  included  in  each  figure  for  comparison,  with  h^ 
being  taken  as  the  average  of  the  depth  at  the  last  grid  point  before  breaking 
and  at  the  first  grid  point  after  breaking.  In  each  case,  the  numerical 
results  provide  an  adequate  representation  of  the  exact  solution. 

4. 5.  Application  to  Offshore  Islands 

Due  to  a  complex  combination  of  wave  breaking,  refraction  and 
diffraction,  the  wave  field  in  the  vicinity  of  an  offshore  island,  either 
natural  or  man-made,  is  extremely  complicated.  For  pure  refraction  Pocinki 
(1950)  has  provided  a  solution,  and  for  nonbreaking  waves  .Jonsson  et  al. 

(1976)  and  Smith  and  Sprinks  (1975)  have  developed  solutions  which  include 
refraction  and  diffraction.  However,  an  island  will  generally  possess  a  surf 
zone  and  the  shadow  zone  created  behind  the  island  will  persist  for  longer 
distances  than  for  the  nonbreaking  wave  case  due  to  loss  of  energy  near  the 
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Downwave  of  the  breakwater,  the  transmission  coefficient  is  given  by 


k-j.  =  |A(x,y)|.  In  Figure  4.9a,b,  several  contour  plots  of  are  provided  for 
breakwaters  with  dimensionless  widths  of  k£  =  5  and  1  Ott  .  Contrary  to  the 
situation  shown  in  Penney  and  Price  (k£  =  10n,  their  Figure  7)  which  neglected 

the  interaction  between  the  diffraction  patterns  from  each  end  of  the 

breakwater,  the  minimum  value  of  kT  does  not  occur  directly  behind  the 
breakwater  but  off  to  each  side  along  a  parabolic  curve  due  to  the  development 

of  a  short  crested  pattern  of  intersecting  waves  in  the  lee  of  the  breakwater 

or  island. 


As  a  non-conservat ive  measure  of  the  extent  of  the  shadow  behind  an 

island,  we  will  use  the  centerline  value  (y  =  0).  Figure  4. 10  shows  k^ 

*  J  nkx-t 

vs.  [  -j .  From  the  figure  the  transmission  coefficient  increases 

asymptotically  to  unity  for  large  x  as  the  waves  diffract  around  the  island. 

2 

(9.3  ki. ) 

As  an  example,  at  a  distance  given  by  kx  =  — — - —  ,  the  wave  height 

recovers  to  90%  of  its  original  value.  Expressing  x  in  terms  of  the  width  of 
the  island,  x  =  Xgg2£,  where  x^q  is  the  dimensionless  distance  to  the  point 
where  90%  of  the  original  amplitude  is  recovered,  we  obtain  x^q  =  [  ]kZ . 

The  distance  xqfl  depends  on  the  relative  size  of  the  island  width  to  a 


wavelength;  the  wider  the  island  the  longer  it  takes  for  the  wave  height  to 
recover.  Note  that  for  kT  =  0.75,  the  distance  is  far  shorter,  x  =  x7s(2i) 


(3.32) 


k£(2?.),  which  is  about  8  times  closer  to  the  island.  In  practice,  a 


factor  of  two  or  three  should  be  applied  to  x,  since,  as  shown,  the  minimum 


value  of  k.j.  is  not  directly  behind  the  island. 


4.7.  Conclusion 

A  method  for  Including  the  effect  of  wave  breaking  in  parabolic  models 
for  combined  refraction  and  diffraction  has  been  developed  and  shown  to  give 
good  agreement  to  laboratory  results  for  wave  heights  acr  the  surface. 


the  model  to  the  other,  instead  of  infinity.  For  simple  geometries, 
analytical  solutions  are  available  as  a  third  way  to  compute  the  far-field 
wave  heights.  This  latter  course  is  chosen  here  to  examine  the  waves  behind 
an  island  approximated  by  an  offshore  breakwater. 

By  superposing  the  solutions  of  the  LaPlace  equation  for  normally 
incident  waves  propagating  past  a  half-infinite  breakwater,  Penney  and  Price 
(1950)  showed  results  of  the  transmission  coefficient  for  an  island  with  a 
width  of  10  wavelengths.  A  similar  superposition  was  done  by  Liu  and  Mei 
(1976)  using  a  similar  parabolic  equation  to  the  one  treated  here.  An  island 
with  a  small  width  to  length  ratio  was  examined  by  Mei  and  Tuck  (1980).  Here 
we  develop  an  exact  solution  of  the  parabolic  approximation  for  the  normally 
incident  wave  field  behind  a  short  breakwater,  in  order  to  approximate  the 
effect  of  an  island  totally  blocking  the  incident  wave.  The  solution  is 
provided  in  Appendix  IT. 

Defining  a  =  /k/nx  (y+H )  and  8  =  /k/irx  (y-2.),  where  k  is  the  wave  number 
of  the  normally  incident  wave,  2JL  as  the  width  of  the  breakwater,  and  x  as  the 
downwave  direction,  the  normalized  wave  amplitude  is 

A(x,y)  =  1  -  (~)  [C(a  )  -  C(8  )  +  i(S(a)  -  S(S  ) )  ]  (4.19) 

where  C(x),  S(x)  are  the  Fresnel  cosine  and  sine  integrals: 
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where  a  was  taken  as  0.2.  This  filter  successfully  eliminated  noise  but 
interacted  with  the  breaking  wave  model  to  the  extent  that  the  overall  form  of 
the  solution  differed.  Instead,  we  employ  a  scheme  in  which  the  model 
searches  for  anomalously  large  values  of  A  occurring  along  grid  rows  with  a 
contracting  surfzone  and  reduces  the  large  values  in  magnitude  to  a  value 
equal  to  the  water  depth.  The  plots  in  Figures  4.6-8  present  unfiltered  data 
in  order  to  give  an  idea  of  the  magnitude  of  the  remaining  noise  and  its 
effect  on  the  overall  contours. 


4.6.  The  Far-Field  of  the  Island 

The  persistence  of  the  shadow  in  the  wave  field  behind  the  island  is 
often  an  important  variable,  particularly  when  siting  an  island  near  developed 
shorelines.  In  order  to  determine  the  wave  field  many  wave  lengths  past  the 
island,  several  procedures  may  be  followed.  First,  the  model  can  be  run  for 
longer  distances  than  shown  in  Figures  4.6-8.  Alternatively,  the  numerical 
model  can  be  used  as  input  to  an  analytical  model,  resulting  in  a  hybrid 
model.  This  can  be  done  using  a  convolution  integral  developed  in  the 
Appendix  to  determine  A(x,y)  over  a  region  of  constant  depth. 

_ _ ,  _i  l  .  i  k(y~p-2- 

A(H,y)  7M7+Te  4J  f(Oe  2(x  °  d4  (4.18) 

where  f(C)  is  the  variation  of  A(x,y)  along  x  =  the  last  grid  row  in  the 
numerical  model.  This  integral  is  limited,  of  course,  to  the  width  of  the 
numerical  model  and  hence  the  Integral  is  taken  from  one  lateral  boundary  of 
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where  Is  the  maximum  height  of  the  conical  island  above  the  bottom  and  XA 
and  Ya  are  the  semi-major  axes  of  the  ellipse  measured  at  the  bottom.  The 
island's  beach  slopes  vary  with  position  around  the  island;  but  on  the  semi¬ 
major  axes,  it  is  h^/X^  and  h^/Y^. 

For  the  computations,  hb  =  20  m  and  X^  and  YA  vary.  The  incident  wave 
has  an  amplitude  of  1  m  and  a  period  of  10  seconds.  The  wave  direction  is 
assumed  to  be  along  the  x-axis  and  the  problem  is  symmetric  about  the  x-axis, 
so  only  half  the  wave  field  is  computed  for  this  case.  Figures  4.6,  4.7  and 
4.8  show  contour  plots  of  the  instantaneous  water  surface  and  the  transmission 
coefficient  around  a  nearly  circular  island  (XA  =  60,  YA  =  80  m),  a  wide 
island  (X^  =  60,  YA  =  160  m)  and  a  long  narrow  island  (XA  =  160,  YA  =  60  m). 
From  the  figures,  which  show  approximately  four  wavelengths,  it  is  clear  that 
the  shape  of  the  island  is  extremely  influential  on  the  wave  field  behind  the 
island,  with  the  island  with  the  smallest  aspect  ratio  creating  the  smaller 
shadow.  Additionally,  the  presence  of  the  island  creates  a  "bow— wave"-like 
disturbance  which  spreads  laterally  in  the  downwave  direction.  The  refraction 
of  the  waves  on  the  front  side  of  the  island  creates  a  region  of  high  wave 
height  on  the  shoulders  of  the  islands  for  all  three  cases.  Further  focusing 
of  the  wave  behind  the  island  occurs  due  to  refraction  by  the  island's 
bathymetry. 

In  practice,  the  shift  from  broken  wave  amplitudes  to  unbroken  wave 
amplitudes  along  a  y-grid  row  Introduces  a  numerical  discontinuity  in  the 
solution.  This  discontinuity  introduces  noise  in  the  computed  solution  which 
then  propagates  into  the  downwave  region.  The  noise  is  manifested  as  a  stable 
Nyquist-wavenumber  oscillation  with  an  amplitude  of  about  one  to  two  percent 
of  the  initial  wave  amplitude.  An  attempt  was  made  to  numerically  filter  the 
oscillation  .luring  model  calculations  by  the  use  of  a  dissipative  interface  in 
the  surfzono  on  the  updated  row: 
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island  boundary.  Several  numerical  (two-dimensional)  examples  of  waves  in  the 
vicinity  of  islands  have  been  computed  here  using  Eqns.  (A.  11)  and  (A.  14)  and 
including  wave  breaking. 

The  inclusion  of  dry  shoreline  boundaries  or  vertical  walls  in  the 
computational  domain  of  a  wave  model  generally  requires  a  modification  of  the 
domain  which  can  include  the  specification  of  stair-step  boundaries  in  the 
case  of  complicated  obstacles.  However,  if  the  surface-piercing  obstacle 
consists  of  a  sloping  structure  which  is  expected  to  be  surrounded  by  a 
surf none,  a  major  simplification  may  be  introduced.  By  replacing  the  surface¬ 
piercing  island  by  a  shoal  with  a  flat  top  and  covered  by  a  thin  layer  with  a 
depth  on  the  order  of  a  centimeter,  the  entire  island  area  may  be  included  in 
the  unmodified  computational  grid.  In  this  "thin  film”  model,  wave  breaking 
then  reduces  wave  height  across  the  surfnone  to  a  small  value  at  the  "real” 
shoreline,  after  which  further  breaking  reduces  the  height  of  a  wave 
propagating  over  the  island  to  approximately  Y  times  the  film  depth.  This 
wave  is  transmitted  beyond  the  island  and  plays  no  role  in  the  downwave  region 
due  to  its  small  height.  This  approach  alleviates  the  need  for  internal  grid 
boundaries  unless  reflective  structures  are  to  be  included. 

The  islands  in  the  present  computations  have  an  elliptical  planform,  and 
are  located  in  water  of  constant  depth  hQ  =  10  m.  The  depth  contours  of  the 
islands  are  given  by 


(A. 16) 
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The  "thin  film"  model  provides  a  convenient  means  to  predict  the  wave 


heights  in  the  vicinity  of  an  island  without  the  computational  necessity  of 
including  internal  boundary  conditions  to  the  model.  Islands  are  treated  as 
extremely  shallow  shoals  and  wave  breaking  insures  that  these  "thin  film" 
shoals  produce  wave  patterns  equivalent  to  those  for  a  surface-piercing 
island. 

The  shadow  thrown  by  islands  in  the  downwave  direction  persists  for 
distances  equivalent  to  many  island  widths,  depending  on  the  dimensionless 
island  width  (k i).  An  analytic  model  for  constant  depth  is  developed  to 
illustrate  the  diffracted  wave  field  behind  the  island. 

There  exists  a  general  absence  of  data  on  the  characteristics  of  breaking 
waves  propagating  either  normally  or  obliquely  over  arbitrarily  uneven 
topography.  Thus,  many  of  the  predictions  of  the  present  model  are 
necessarily  qualitative  in  nature,  and  require  further  verification  which 
depends  on  the  future  availability  of  data. 


Appendix  4.1.  Analytic  Solution  for  the  Constant-Depth  Half-space 

Over  a  flat  bottom,  neglecting  dissipation,  the  parabolic  equation  for 
t!ie  wave  amplitude,  A(x,y),  behind  the  island  is 


2ik 


0 


'  is  the  wave  number.  This  equation  may  be  used  to  extend  the  results 
•r  / i ous ly  described  finite-difference  calculations  to  the  farfield  of 
>r  linilir  obstacle).  The  boundary  condition  (obtained  at  the 
•  r(ir  eunputat ional  grid)  may  be  taken  as 


A(0,y)  =  f (y)  ;  |y|<“ 

where  now  x  =  0  signifies  the  start  of  the  half-space  to  be  considered.  Using 
the  Fourier  transform  pair: 

CO  oo 

A  =  /  A(x,y)e  lXydy  ;  A  =  j  A(x,X )e1XydX  , 

—CO  — oo 

we  obtain  an  equation  in  x  alone, 

3  A  .2; 

2ik  - —  -  X  A  =  0 
3  x 

where 

oo 

A(0,X)  =  f(X)  =  ■  f,v)e  iAydv 

Solving,  we  obtain 

w  ,  ^  ,  >.  X  "x/2ik 

A(x,X  )  =  A(0,X  )e 

This  solution  may  be  used  in  either  continuous  form  or  in  a  discretized 
form  using  FFT  methods  to  extend  the  finite-difference  calculations.  Here,  we 
consider  a  special  case  of  a  breakwater  of  length  21  situated  on  the  y-axis  at 
x  =  0.  For  the  general  case  of  waves  obliquely  incident  on  x  =  J,  f(y)  may  be 
written  as 


J  y  }<«- 

I  y  I  >a 
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where  XQ  is  the  y-component  of  the  wavenumber  vector 


XQ  =  ksin8 


Transforming  f(y)  gives 


„  2sin(X  -X )l 

A(0,X)  =  2tt  6  (Xq-X  )  - 

o 


where  <5  is  the  Dirac  delta  function. 


Taking  the  inverse  transform  of  A  gives 


2sin(X  -\)l  2 


.  /  \  1  f  r .  '  o  i  -iX  x/2k  iXy 

A(x,y)  =—  J  [2tt6(X  -X)  -  — -  . - J  e  e  *dX 


(X  -X) 

O 


-iX  x/2 k  iX  y  i— — r 

o  oJ  /  k  -m/4  t 

=  e  e  J—  e  J 


/*  elA°Celk(y‘4)2/2’‘dc 


where  the  first  term  represents  an  undisturbed,  obliquely-oriented  plane  wave 
and  the  second  term  contains  the  disturbance  due  to  the  breakwater.  Setting 
XQ  equal  to  zero  (normal  incidence)  and  introducing  the  definitions  of  the 
Fresnel  cosine  and  sine  integrals  yields  equation  (4.17).  The  solutions 
obtained  using  the  parabolic  approximation  are  similar  in  form  to  the  form  of 
the  Fresnel  approximation  in  optics. 
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